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Abstract 

The variance of the local density of the pair contact process with diffusion (PCPD) is investigated 
in a bosonic description. At the critical point of the absorbing phase transition (where the average 
particle number remains constant) it is shown that for lattice dimension d > 2 the variance exhibits 
a phase transition: For high enough diffusion constants, it asymptotically approaches a finite value, 
while for low diffusion constants the variance diverges exponentially in time. This behavior appears 
also in the density correlation function, implying that the correlation time is negative. Yet one has 
dynamical scaling with a dynamical exponent calculated to be z = 2. 



I. INTRODUCTION 



An prototypical example for critical phenomena in nonequilibrium statistical physics is 
the absorbing phase transition. This is a transition from an active fluctuating phase with 
a finite particle density to an absorbing state where any dynamics is suppressed. One has 
found rather robust universality classes, e.g. the class of directed percolation (DP) and the 
parity conserving universality class (PC). A member of the DP-class is the pair contact 
process where two neighboring particles may create an offspring on a third lattice site or 
may annihilate each other. 

This model extended by particle diffusion - the pair contact process with diffusion 
(PCPD) - has attracted much interest because it is not known to which universality class 
it belongs. Several possibilities have been discussed: It was found that some exponents are 
very close to those of the PC class more recent investigations however give hints for a 
DP behavior Q]. It was also suggested that the critical behavior of the PCPD defines a 
new universality class , or may depend on the diffusion constant [5]. Analytical results 
are rare in this field and one has to revert to numerical methods like the density matrix 
renormalization group (DMRG) or Monte Carlo simulations. 

This is different for the bosonic description of the model where the exclusion interaction 
- which constraints the number of particles at one site to at most one - is dropped. In 
this case a field theoretic approach due to Howard and Tauber [f| is available. A drawback 
of this approach is that it is not suitable for deciding the universality class of the model 
with particle number restriction. In this paper we show by an exact treatment of the model 
that the diffusion constant and the lattice dimension have considerable impact on the phase 
transition and correlations of the bosonic PCPD. Although the particle exclusion interaction 
is crucial for the behavior of the system this investigation gives some insight on the role of 
diffusion in the PCPD. 

II. MODEL 

We define the following process: On a infinite <i-dimensional cubic lattice particles ('^4') 
are diffusing with rate D, in each spatial direction. Additionally they branch and annihilate: 
k > 1 particles A are created with rate fi out of any set of m > 1 particles (m fixed), and 
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I > 1 particles are annihilated with rate A out of any set of p > I particles (7 fixed): 



mA — > [m + k)A 
pA ^ {p-l)A 
A- & -A. 



(1) 



The number of particles on each lattice site is not restricted - the creation and annihilation 
processes take place on one lattice site. Thus the bosonic representation of the process is 
used. We try to keep the description as general as possible, but as we will see, analytical 
results are available only for few cases. In this paper we investigate the two cases where 
p = m = lorp = m = 2 and arbitrary k and / < p. One special case is the PCPD, where 
m = p = I = 2 and k = 1. 

Following the notation and formalism introduced in [?], O] we define the site occupation 
numbers as n = {n(x)}. Then the time dependent probability vector describing the system 
can be expressed as 

\F(t)) = J2P(n,t)\n) (2) 

n(x) 

where the \n) are the basis vectors spanning the state space and P(n, t) is the probabil- 
ity distribution of the site occupation numbers. The master equation describing the time 
evolution of the probability distribution can then be written as 

d 



Of 



\F(t)) = -n\F(t)) 



(3) 



7i is the stochastic generator of the system, often called as "hamiltonian" due to the analogy 
of the master equation to the Schrodinger equation (in imaginary time) |3|. Let a(x) and 
o(x)' be the space dependent annihilation and creation operators and ra(x) = cr(x)a(x) the 
particle number operator, then the hamiltonian is given by 

d 

H = ~D^2Y^ [a(x)a f (x + k) + a f (x)a(x + k) - 2n(x) 



k=l x 



X 



(at(x)) (p -° (a(x)) p -;QHx)-z + l) 



8=1 



(at(x)) (m+fc) (a(x))--n 



nix 



i=i 



(4) 
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where k = k(&) = (...,0,1,0,.. .) T is the k-th unit space vector. The time evolution of an 
operator b(y) is calculated by 

^b(y) = [H,b(y)}- (5) 
Using the commutator rule [a(x),a^(y)j = 5 x . y we get after straightforward calculations 

!<a(x)> = J D^{(a(x-k)) + (a(x + k))-2(a(x))} (6) 



k=l 

-A/(a(x) p ) +/ifc(a(x) m ) 

d 



^(a(x)a(y)) =^ D ^ { (a(x)a(y — k)) + (a(x)a(y + k)) + (7) 

k=l 

(a(x - k)a(y)) + (a(x + k)a(y)) - 4(a(x)a(y)) } 
-XI { (a(x)a(yf) + (o(x)'a(y)> } 
+/iA;{ (a(x)a(y)" 1 ) + (a(x) m a(y)> } 

((a(x)) 2 ) = 2 J D^{(a(x)a(x-k)> + (a(x)a(x + k)) -2(a(x) 2 )} (8) 
fc=i 

+A/ { (1 + Z - 2p) (a(x) p ) - 2(a(x) p+1 > } 

m+l\ 



<9f 



-/ifc { (1 - k - 2m)(o(x) m ) - 2(a(x) m+1 ) } 

Using («-(x)) = (a(x)) and (n(x) 2 ) = (a(x) 2 ) + (a(x)) this set of coupled difference- 
differential equation allows for the analytical calculation of the time-dependent expectation 
value of the particle density and its autocorrelation in some special cases. 

We restrict to the case p = m where the creation and annihilation processes are balanced 
and an absorbing phase transition can be found. For XI > fik the particles die out expo- 
nentially (p = m — 1) or according to a power law (p = m > 1), while for XI < fik the 
particle density diverges. Here a crucial difference between the description with and with- 
out particle number restriction can be seen: While in the models with exclusion interaction 
the absorbing phase transition is of second order, the bosonic model exhibits a first order 
transition. 

In analogy to the exclusion model we call the rate which divides the two different behaviors 
the "critical" rate, which from Eq. © can be read off as 

A c = fik/l (9) 

for given /i. For this rate the particle density is constant for all times (a(x, t)) = po (f° r 
homogeneous initial conditions), as can be seen from Eq. © which reduces to a diffusion 



equation. Thus the interesting quantity is the variance a 2 = (n(x) 2 ) — (n(x)) 2 which we 
shall investigate in what follows. 

Eliminating p and A in Eqs. ©-(JEJ) one gets 

|(a(x)> =D «a(x - k)> + (a(x + k)> - 2(a(x))} 

k=l 

d d 

— (a(x)a(y)) = D { (a(x)a(y - k)> + (a(x)a(y + k)) + 



(a(x - k)a(y)) + (a(x + k)a(y)) - 4(a(x)a(y)) } 
|((a(x)) 2 ) =2D £ { (a(x)a(x - k)> + (a(x)a(x + k)> - 2(a(x) 2 )} 



(10) 



fc=i 

+ lik(k + l)(a(x.) m ) 

We see that this set of equations is only closed for the cases m = 1 or m = 2. 
In the case of a vanishing diffusion constant, D = 0, the lattice sites are independent of 
each other. Thus the description of the process reduces to the zero dimensional case d = 0, 

!:<a( x )> =0 



dt 

((a(x))*> =^(A; + /)(a(x) m ), 



dt 

and has to be treated separately. 

A. Contact process with diffusion, m = 1 

Here, only I = 1 is possible. Additionally by rescaling /i we may fix = 1. This case has 
already been considered in as a model for clustering of biological organisms Q]. For 
convenience we summarize the main results here. 

For D = or d = Eq. (jllj) directly yields (a(x) 2 ) = cq + ci £ and thus the variance 
diverges. For D ^ the fluctuations of the particle density diverges for dimensions d < 2 
while they remain finite for d > 2, 

ci r d / 2+1 d < 2 
(a(x) 2 ) = <J C2 hit d = 2 (12) 
c 3 + c 4 r d/2+1 d > 2 
where t > 1 and cq, . . . , C4 are positive constants. 



B. Pair contact process with diffusion, m = 2 



We now derive analytically the late time behavior of the solution for m = 2. 
For D = or d = Eq. fllTJ) yields 

(a(x 2 )) = p 2 exp(t/T), 

1 (13) 
pk(k + I) 

The variance diverges exponentially in time as opposed to m — 1 where the divergence is 
linear. Only for times small compared to r the variance Eq. IjlHj) grows linearly. 

For D ^ we get the solution by applying Fourier- and Laplace-transformations. We 
also present the crossover from short to late time behavior, which has to be calculated 
numerically. 

First we rescale time by 

t - 4 (14) 

and define 

F x (r, t) =(o(x)o(x + r)) = (n(x)n(x + r)) - <J P , (n(x)) 

uk(k + l) ^ 

a = . 

2D 

The parameter a is a measure for the weighting of reaction rates to diffusion, small a 
corresponds to dominant diffusion, while large a corresponds to dominating reaction rates. 
In what follows we consider only translational invariant inital conditions, in which case 
F x (r,t) is independent of x. Using Eq. (|10p we get the following difference-differential 
equation for F: 

d 

Fir -n = ' 



9 -F(r, t)=J2 i F (r ~ M) + F(r + k, t) - 2F(r, t)} + 5 r , aF(0, t) 



k=i 



(16) 



J2 A kF(r,t) + 5 r , o aF(0,t) 



k=l 

'2 



where is the discrete Laplacian concerning the fc-th component. The variance o is 
related to F as follows 

a(t) 2 = F(0,t)+p -p 2 . (17) 

Here, we see that there is no qualitative difference between parity conserving models (k 
and / even) and non-parity conserving models — models with different k and I differ only 
by different creation and annihilation rates. 



This kind of equation can be solved using the Fourier-transformation: 

/(q,t) = ^e-^F(r,t), F(r,t) = J -^ e *«7(q,f). (18) 

We get 

^/(q, t) = -w(q)/(q, t) + aF(0, t), (19) 
with the dispersion relation w(q) = —2 Xlfc=i ( cos ilk) ~ 1)- Integration yields 

/(q, f) = e-^ |/(q, 0) + a jf drF(0, r)e w ^| . (20) 

As initial condition we choose a Poisson-distribution F(r, 0) = pi so that /(q, 0) = <5 qi oPo- 
Thus we get 

F(r,t)=p2 + a / drF(0,r)6(r,t-r) (21) 
Jo 



with 

&(r,t) = [ J^ e -™(<i)t+w 

V > J W (22) 
=e- 2dt I ri (2t)-...-I rd (2t) 

where I r (t) is the modified Bessel function of order r. The dimension d is now just a 
parameter which can formally take real values. Although this is not physical it allows for 
the investigation of the dependence on the dimension. 

For r = the long-time behavior of the solution of the Volterra integral-equation Eq. (|21|) 
with the function b(t) given by Eq. (}2*2*j) is known from the mean spherical model |l5j. In this 
context a plays the role of the temperature. This analogy enables us to use known results 
rom the spherical model. Eq. (|2*T|) can be solved using temporal Laplace transformation 

3, 

/■oo 

F(p) = / dte- pt F(0,t). (23) 
Jo 

We get 

p(p) = d + a F(p)b(p) 
V 

* HP) = n Pl ~ h(w (24) 
p[l — ab{p)) 
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for late times -F(0, t) is given by the behavior of b(p) for small p, which crucially depends 
on the dimension d (see for example jl3^: 



b(p) 



Ai 



(4n)-*/*r(l-d/2)p-V-- d M d<2 

2A X - (47r)- d / 2 |r(l - d/2)\p d l 2 - 1 2<d<4 

2At - AA 2 p d>4 
d d q 1 



(25) 



(2Tv) d (2w(q)) k 

This results in different behavior of -F(0, t) as we shall see in the next sections. 

For all even integral dimensions d — 2, 4, . . . logarithmic corrections arise whose investigation 

goes beyond the scope of this paper. 



1. d < 2 

As for d < 2 the quantity b(p) diverges for p — > the denominator of Eq. (|24jl has always 
a zero for p ^ 0, so that F(p) has a pole at a positive value p — 1/r. A pole of the laplace 
transform corresponds to exponential behavior of the original function and we get 

F(0,t) oc e i/r . (26) 



t— >oo 



For d = 1 the exact expression of b is known 



~ m = 7m (27) 

which yields 

r d =i = ; 1 • (28) 

For any finite value of a the time scale r is finite but diverges if a \ 0. This is in analogy 
to the spherical model, where in one dimension the critical temperature is zero. 



In order to investigate how the predicted asymptotic behavior for large times is ap- 
proached we have performed a numerical integration of -F(0, t), shown in Fig. ^ For details 
of the numerical calculation see ^(|, where a similar integral equation is calculated. We 
see that the asymptotic behavior is approached quickly and the solution Eq. (|2*S|) is a good 
approximation for times t > 1. 
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2. 2 < d < 4 



For d > 2 the quantity b(p) shows qualitatively different behavior: It approaches the 
finite value 2A 1 for p — > 0. Therefore the F(p) has a pole for positive p only for ct larger 
than a critical value given by 

a. = ^, (29) 

which is identical to the critical temperature in the spherical model. Thus we find a phase 
transition in the behavior of the autocorrelation F(0, t): For a > a c (low diffusion constant) 
we recover the exponential divergence 

F(0,t) oc e t/T . (30) 

with a time scale 



— (31) 
c 2 a ) 

with the reduced control parameter 

a' = ^— ^, (32) 

and c 2 = (47r) _d / 2 |r(l — d/2)\. This time scale diverges if we approach a \ a c . 

For a < a c (high diffusion constant) the pole of F(p) vanishes and F(0, t) asymptotically 
approaches a finite value 

F^ = lim F(0, t) = p2 °, > pi (33) 
t^oo i — a/a c 

which diverges if we approach a S a c . 

Therefore a suitable order parameter for this phase transition is F^ 1 which decreases 
linearly to zero for a S a c and is equal to zero for a > a c . 

For a = a c we get 

F{P) ~ \F(l-d/2)\a c p^ (34) 

which results in a power law 

F(0,t) (xt d/2 - 1 (35) 
and hence in a power law divergence of the variance. 
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3. d > 4 



For d>4we find qualitatively the same behavior as for 2 < d < 4. Like in the previous 
case, F(p) has a pole at a positive p only for values a > a c = l/(2Ai). For a > a c the time 
scale of the exponential increase is given by 

a' 



Mi*,' ■ m 

The difference to the case 2 < d < 4 is, that this time scale is now independent of the 
dimension d, indicating that we are in the mean field region. 

For a < a c , F(0,t) approaches the asymptotic value given by Eq. (JSHJ- 

For a = a c we get 

*® ^ W (37) 

which results in a power law 

F(0,t)oct (38) 
These results are summed up in the phase diagram Fig. |3J 



III. SPATIAL CORRELATIONS 



In the mean field regime (d > 4) the behavior of the correlation function F(r, t) can be 
calculated analytically in the limit of large r and t. As derived in the appendix we get: 



P0 a 
47T d / 2 ]0 



„2-d ■ 



r li-ix 



' it 



F{T,t)-p 2 ={ 



Pi 



64A 2 7r d / 2 



r 4-d ^ d 



i) 



a' < 



a' = 



(39) 



i (8 ,)(^>/^- 2 >/3 (i) {d ' 4)/2 (l) (1 "" )/2 ex P (* A - r/0 i » ct > 



where a' = (a — a c )/a c is the reduced control parameter, T is the incomplete Gamma 
function and \1/ is a scaling function defined by 

r(f-M) 



ty(d : u) 



dz- 



(40) 



Above the critical point the correlations diverge; the time scale r is given by Eq. (j3*U|) and 
the correlation length by 

^ (41) 



Q 
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Interestingly, as for a c > the correlations increase with time, what in usual dynamical 
critical phenomena would be called the correlation time is negative while the correlation 
length is positive. For a' < the dependence on r 2 /t directly shows that the dynamical 
exponent is z = 2. For a' > the time scale r is the square of the length scale £, therefore 
also in this case the dynamical exponent is z = 2. 

No analytical solution is available in the case 2 < d < 4, thus we evaluate the integral 
(12 lj) numerically Fig. 0] shows the spatial dependence of the correlation function along the 
axis r = (r, 0, ...), for a < a c ; Fig. |5] shows the case a > a c . A collapse of the calculated 
data points is achieved if we assume the following functional dependence: 



where f\ is a scaling function and f2{d,r) is a function that only depends on d and r. This 
result is in qualitative agreement with the previously derived formula for d > 4. 

IV. DISCUSSION 

Apart from the spherical model this phase transition is related to a much simpler model: 
On a rf-dimensional cubic lattice non interacting particles are diffusing with rate D and 
additionally at site x = particles may branch A — > 2 A with rate a' = aD. The equation 
for the time evolution of the particle density (n(x, t)) is just given by Eq. ()16|). We can adopt 
the solutions for F(r,t) by substituting the initial condition by — > p . In particular we 
recover a phase transition for the particle density at the origin. While in the original process 
it is rather complicated to understand the physical meaning of the behavior of the second 
moment, in this model we understand the behavior of the first moment: For d = 1 diffusion 
does not suffice to spread the particles on the lattice fast enough and the particle density at 
x = diverges for any given parameters. For higher dimensions additional spatial directions 
are accessible to spread particles and as a consequence the particle density at x = remains 
finite for high enough diffusion constant D. 

The fact that the autocorrelation function is diverging while the particle density remains 
constant allows some conclusions concerning the distribution function for the particles p(n) 
for late times. On the one hand, if (n) = ^2 n np(n) is finite then for large n the distribution 




(42) 
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function p(n) < c\rT^ with (3 > 2. On the other hand, if (n 2 ) = 2~2 n n2 P( n ) ^ s infinite 
then for large n the distribution function p(n) > c 2 n _/3 with (3 < 3 with some positive 
constants ci, c 2 . Thus the distribution function follows for large n a power law p{n) oc n^ 13 
with 2 < /3 < 3. 

In summary, we have shown that for d > 2 the bosonic PCPD exhibits a phase transition 
for (a(x) 2 ) and thus for the autocorrelation function cr(t) 2 = (n(x) 2 ) — (n(x)) 2 = (a(x) 2 ) + 
(n(x)) — (n(x)) 2 . The order parameter F^ 1 decreases linearly to zero for a / a c and is 
equal zero for a > a c , where a is proportional to the ratio of the reaction rates and the 
diffusion constant. Thus diffusion has big influence in this process, it must be high enough 
in order to avoid a divergence of the autocorrelation. 

We have also shown that the critical properties of this process are related to the mean 
spherical model. As the spherical model is a model for magnetism this analogy is rather 
intriguing and the question arises whether it is just accidental. 



For the meanfield case d > 4 the solution of F(r,t) in the limit of large r and t can 
derived analytically, as presented in what follows. 
With the definition 



Appendix 



G(r,t) = F(v,t) - pi 



(A.43) 



the integral equation can be transformed to 




(A.44) 



Using a Laplace transformation we get 



G(r,p) 



aG(0,p)b(r,p) + apl 



2 b(r,p) 



(A.45) 



P 



setting r = determines G(0, p) which yields 




b(r,p) 



(A.46) 



The Fourier transform of this equation is 




1 



1 



(A.47) 
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For the meanfield case, d > 4, 6(0, p) takes the simple form 

~b(0,p) = l/a c - P1 /a, (A.48) 

and we get 

Pc\0i 1 1 , . , , 

»(q,P =- 7 77^ - 7 v A.49 

Here we defined the reduced control parameter a' = (a — a c )/a c and 7 = 4A 2 a. Although 
the function b does not depend on dimension for d > 4, generally a dependence of the 
solution G(r, t) on dimension is still possible as the inverse fourier transform depends on d, 
which does not affect the critical exponents. Using an expansion into partial fractions we 
get 

„, s 2 / 1 1 1 1 1 1 

</(q,P) = Po a --7-7-T- + -TT-yZV—; 77777 1777 + 



a;'w(q)p a' (io(q) + a'/ 7) p — ac'/'y 7w(q) (iw(q) + ck'/t) P + w ( c l) 

(A.50) 

The inverse Laplace transform of this expression reads 

^ t] = Pla (-^) + ^KqTT^) ^ ( 7 ') + 7^(q) Md) + aV7) ^ 

(A.51) 

Although the second term is not Laplace transformable for a' > this result is correct and 
can be derived by transforming the function H(r,t) = exp + e)t) G{r,t) with e > 0. 

The inverse fourier transform of the first term of (|A.51|) is 

f d d q iqr 1 [°°. f d\ ■ 

/ 7 — he —r^= / dx / - — -f-r exp (-w(q)a;) e q 

y (27r) d w (q) y y (2^ Pl KH) > 

POO 

= / dxe- 2dx I ri (2x) ■ ... -I rd (2x) 
Jo 

J dx (47rx) _d/2 exp (^--^j 



(A.52) 



|r|>l 

_ r(|-i) 



4 7r d/2 

For the long time limit the second term contributes significantly only for a' > 0. Defining 
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b 2 = C//7 we get for this case: 

f d " q - = [°° dx I exp (- (w(q) + b 2 ) x) e 



f 

Jo 



dxexp(-b 2 x)e- 2dx I ri (2x) I rd (2x) 



6 2 <l,|r|>l 



(4tt) - d/2 J™ dx x~ d/2 exp (- ^ - 6 2 x 



(47T) 



-d/2 i_ 



/ 2\l~d/2 /-oo / l2 2 

f^J J dzz d ^ 2 exp U - 6 ' 



4z , 

2\ i-d/2 (A.53) 



(2 \ i— a/ z 
(br)'" 2 ' 1 K m _, (6r) 

i ^'V" ,/4 r(1 - J)/2ex /_ R 



2 (cH-l)/2 7r (d-l)/2 



7/ V V 7 



where K d / 2 -\ is the modified Bessel function of second kind. 
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For a' > the third term is transformed to: 

d d q e iqr exp(-w(q)t) 
(27i) d w(q) (w(q) + b 2 ) 

Ha f d ^ ( f f , e^exp(- W (q)t) 

/■°° r 00 /" d d q 

= J dxexp(-b 2 x) J dy J -^—^ exp (-w(q)y) e iqr 

/•CO /»oo 

= / dxexp(-6 2 x) / dye- 2 *'/ ri (2y) • . . . • I r<f (2y) 

JO Jt+x 

roo poo / 2 \ 

« Un)- d/2 dxexp(-b 2 x) dyy- d/2 exp( 

M>i Jo ' Jt+x V A VJ 

= (4 7 r)- d / 2 t 1 - d / 2 ^ dx exp(-6 2 a;) dz exp( . 



5i4^y ^^(-^^ d^^exp(-z)-^-j exp(--)- 

r 2-d roo / / ^ \ / ^ r 2\ / r 2\ a! / 2 - 1 r 2 



dz z d/2 - 


_1 exp(- 








V 


















v4t 



exp(- 



r 2 ~ d ,_ 2 („fd _\ „ / d r 2 \ 6 -2 / r 2 \ d//2 ~ ] „ r 2 



(A.54) 



- 4,^ { T {2- l )- T {-2- 1 'Tt)--{T t ) 

with the incomplete Gamma-function T(a,x) = dtt a_1 e _i . 

For a' < the same result can be derived by substituting x — > — x. 
At the critical rate a' = equation (|A.49|) reduces to 

~9(^P) = ^\ * v (A.55) 
The inverse Laplace transform of this expression is given by 



C- 1 ( -, 1 . .. ,t) = f dr [dT'C^l 1 —-,t' 

= dr dr' exp(— lu(q)r'). 
Jo </o 

The necessary conditions for this equality are fulfilled [l^j : 

lim (^' Pt J drexp(— w(q)r)^ =0 
hm ( e_P *^ dr ^ dr / exp(-w(q)r')^ =0. 
15 



(A.56) 



(A.57) 



This yields 

-,2^, ft ft r A d 



G(M)= T/ dT I dT 7(^ exp( - u,(q)r,) 

f dr f dr' e- 2dx J ri (2x) • . . . • I rd {2x) 
Jo Jo 
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« ^ fdr rdr / (47rr')- d/2 expf-^ 
|r|»i 1 Jo Jo V 4 W 

n 2^ /„2\-d/2+l ,t ,00 (A.58) 



( 2\ /-t /-oo 

4t 



7(47r) J / 2 V 4 



>2 



4t 

Thus we get the expressions Eq. ()39|) in the limit of large r and t. 
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FIG. 1: Numerical calculation of F(0,t) for d = l,a = 2,po = 0.1. The dashed line shows the 
theoretical predicted slope 1.2071. 




FIG. 2: Numerical calculation of F(0, t) for d = 3, po = 1, a = 2 < a c (solid line) and a = 4.2 > a c 
(dashed line). The dotted line shows the theoretical predicted slope r w 30.4, the dashed-dotted 
line the theoretical predicted asymptotic value F(0, t = oo) ~ 2.02. 
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FIG. 3: The phase diagram of the system for fixed diffusion constat D: In the limit of t — > oo 
for A < fik/l the stationary density p* is zero while it diverges for A > pk/l. For A = pk/l the 
density is constant, p* = po, and the variance function is bounded for a < a c , while it diverges 
exponentially for a > a c and algebraically for a = a c , where a = pk(k + 1)/ (2D). 




FIG. 4: Numerical calculation of F(r = (r, 0, ...), t) for d = 3, a = 2 < a c , Po = 1 and times 
t = 50,100,150,... ,400. 
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FIG. 5: Numerical calculation of F(r = (r, 0, ••■),£) for d = 3, a = 6 > a c , po = 1 and times 
t = 50, 100, 150, ...,400. 
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